library(tidyverse)
library(Seurat)
library(cowplot)
library(ComplexHeatmap)
library(circlize)
library(GeneOverlap)
library(gprofiler2)
library(ggrepel)
library(ggplot2)
library(muscat)
library(purrr)
library(limma)
library(scran)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ── ✔ dplyr 1.1.0 ✔ readr 2.1.4 ✔ forcats 1.0.0 ✔ stringr 1.5.0 ✔ ggplot2 3.4.1 ✔ tibble 3.1.8 ✔ lubridate 1.9.2 ✔ tidyr 1.3.0 ✔ purrr 1.0.1 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ── ✖ dplyr::filter() masks stats::filter() ✖ dplyr::lag() masks stats::lag() ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors Attaching SeuratObject Attaching package: ‘cowplot’ The following object is masked from ‘package:lubridate’: stamp Loading required package: grid ======================================== ComplexHeatmap version 2.14.0 Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/ Github page: https://github.com/jokergoo/ComplexHeatmap Documentation: http://jokergoo.github.io/ComplexHeatmap-reference If you use it in published research, please cite either one: - Gu, Z. Complex Heatmap Visualization. iMeta 2022. - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 2016. The new InteractiveComplexHeatmap package can directly export static complex heatmaps into an interactive Shiny app with zero effort. Have a try! This message can be suppressed by: suppressPackageStartupMessages(library(ComplexHeatmap)) ======================================== ======================================== circlize version 0.4.15 CRAN page: https://cran.r-project.org/package=circlize Github page: https://github.com/jokergoo/circlize Documentation: https://jokergoo.github.io/circlize_book/book/ If you use it in published research, please cite: Gu, Z. circlize implements and enhances circular visualization in R. Bioinformatics 2014. This message can be suppressed by: suppressPackageStartupMessages(library(circlize)) ======================================== Loading required package: SingleCellExperiment Loading required package: SummarizedExperiment Loading required package: MatrixGenerics Loading required package: matrixStats Attaching package: ‘matrixStats’ The following object is masked from ‘package:dplyr’: count Attaching package: ‘MatrixGenerics’ The following objects are masked from ‘package:matrixStats’: colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins, rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans, rowWeightedMedians, rowWeightedSds, rowWeightedVars Loading required package: GenomicRanges Loading required package: stats4 Loading required package: BiocGenerics Attaching package: ‘BiocGenerics’ The following object is masked from ‘package:limma’: plotMA The following objects are masked from ‘package:lubridate’: intersect, setdiff, union The following objects are masked from ‘package:dplyr’: combine, intersect, setdiff, union The following objects are masked from ‘package:stats’: IQR, mad, sd, var, xtabs The following objects are masked from ‘package:base’: anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply, union, unique, unsplit, which.max, which.min Loading required package: S4Vectors Attaching package: ‘S4Vectors’ The following objects are masked from ‘package:lubridate’: second, second<- The following objects are masked from ‘package:dplyr’: first, rename The following object is masked from ‘package:tidyr’: expand The following objects are masked from ‘package:base’: expand.grid, I, unname Loading required package: IRanges Attaching package: ‘IRanges’ The following object is masked from ‘package:lubridate’: %within% The following objects are masked from ‘package:dplyr’: collapse, desc, slice The following object is masked from ‘package:purrr’: reduce Loading required package: GenomeInfoDb Loading required package: Biobase Welcome to Bioconductor Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'. Attaching package: ‘Biobase’ The following object is masked from ‘package:MatrixGenerics’: rowMedians The following objects are masked from ‘package:matrixStats’: anyMissing, rowMedians Attaching package: ‘SummarizedExperiment’ The following object is masked from ‘package:SeuratObject’: Assays The following object is masked from ‘package:Seurat’: Assays Loading required package: scuttle
library(future)
#for 200gb ram
options(future.globals.maxSize = 200000 * 1024^2)
sessionInfo()
R version 4.2.2 (2022-10-31) Platform: x86_64-conda-linux-gnu (64-bit) Running under: AlmaLinux 9.3 (Shamrock Pampas Cat) Matrix products: default BLAS/LAPACK: /hpc/group/pbenfeylab/tmn23/miniconda3/envs/muscat/lib/libopenblasp-r0.3.21.so locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=en_US.UTF-8 LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats4 grid stats graphics grDevices utils datasets [8] methods base other attached packages: [1] future_1.31.0 scran_1.26.0 [3] scuttle_1.8.0 SingleCellExperiment_1.20.0 [5] SummarizedExperiment_1.28.0 Biobase_2.58.0 [7] GenomicRanges_1.50.0 GenomeInfoDb_1.34.8 [9] IRanges_2.32.0 S4Vectors_0.36.0 [11] BiocGenerics_0.44.0 MatrixGenerics_1.10.0 [13] matrixStats_0.63.0 limma_3.54.0 [15] muscat_1.12.0 ggrepel_0.9.3 [17] gprofiler2_0.2.1 GeneOverlap_1.34.0 [19] circlize_0.4.15 ComplexHeatmap_2.14.0 [21] cowplot_1.1.1 SeuratObject_4.1.3 [23] Seurat_4.3.0 lubridate_1.9.2 [25] forcats_1.0.0 stringr_1.5.0 [27] dplyr_1.1.0 purrr_1.0.1 [29] readr_2.1.4 tidyr_1.3.0 [31] tibble_3.1.8 ggplot2_3.4.1 [33] tidyverse_2.0.0 loaded via a namespace (and not attached): [1] pbdZMQ_0.3-9 scattermore_0.8 [3] bit64_4.0.5 irlba_2.3.5.1 [5] DelayedArray_0.24.0 data.table_1.14.8 [7] KEGGREST_1.38.0 RCurl_1.98-1.10 [9] doParallel_1.0.17 generics_0.1.3 [11] ScaledMatrix_1.6.0 RhpcBLASctl_0.23-42 [13] RSQLite_2.2.20 RANN_2.6.1 [15] bit_4.0.5 tzdb_0.3.0 [17] spatstat.data_3.0-0 httpuv_1.6.9 [19] viridis_0.6.2 hms_1.1.2 [21] evaluate_0.20 promises_1.2.0.1 [23] fansi_1.0.4 progress_1.2.2 [25] caTools_1.18.2 igraph_1.3.5 [27] DBI_1.1.3 geneplotter_1.76.0 [29] htmlwidgets_1.6.1 spatstat.geom_3.0-6 [31] ellipsis_0.3.2 backports_1.4.1 [33] annotate_1.76.0 aod_1.3.2 [35] deldir_1.0-6 sparseMatrixStats_1.10.0 [37] vctrs_0.5.2 ROCR_1.0-11 [39] abind_1.4-5 cachem_1.0.6 [41] withr_2.5.0 progressr_0.13.0 [43] sctransform_0.3.5 prettyunits_1.1.1 [45] goftest_1.2-3 cluster_2.1.4 [47] IRdisplay_1.1 lazyeval_0.2.2 [49] crayon_1.5.2 genefilter_1.80.0 [51] spatstat.explore_3.0-6 edgeR_3.40.0 [53] pkgconfig_2.0.3 nlme_3.1-162 [55] vipor_0.4.5 blme_1.0-5 [57] rlang_1.0.6 globals_0.16.2 [59] lifecycle_1.0.3 miniUI_0.1.1.1 [61] rsvd_1.0.5 polyclip_1.10-4 [63] lmtest_0.9-40 Matrix_1.5-3 [65] IRkernel_1.3.2 boot_1.3-28.1 [67] zoo_1.8-11 base64enc_0.1-3 [69] beeswarm_0.4.0 ggridges_0.5.4 [71] GlobalOptions_0.1.2 png_0.1-8 [73] viridisLite_0.4.1 rjson_0.2.21 [75] bitops_1.0-7 KernSmooth_2.23-20 [77] Biostrings_2.66.0 blob_1.2.3 [79] DelayedMatrixStats_1.20.0 shape_1.4.6 [81] parallelly_1.34.0 spatstat.random_3.1-3 [83] beachmat_2.14.0 scales_1.2.1 [85] memoise_2.0.1 magrittr_2.0.3 [87] plyr_1.8.8 ica_1.0-3 [89] gplots_3.1.3 zlibbioc_1.44.0 [91] compiler_4.2.2 dqrng_0.3.0 [93] RColorBrewer_1.1-3 clue_0.3-64 [95] lme4_1.1-31 DESeq2_1.38.0 [97] fitdistrplus_1.1-8 cli_3.6.0 [99] XVector_0.38.0 lmerTest_3.1-3 [101] listenv_0.9.0 patchwork_1.1.2 [103] pbapply_1.7-0 TMB_1.9.2 [105] MASS_7.3-58.2 tidyselect_1.2.0 [107] stringi_1.7.12 BiocSingular_1.14.0 [109] locfit_1.5-9.7 tools_4.2.2 [111] timechange_0.2.0 future.apply_1.10.0 [113] parallel_4.2.2 uuid_1.1-0 [115] bluster_1.8.0 foreach_1.5.2 [117] metapod_1.6.0 gridExtra_2.3 [119] Rtsne_0.16 digest_0.6.31 [121] shiny_1.7.4 Rcpp_1.0.10 [123] broom_1.0.3 later_1.3.0 [125] RcppAnnoy_0.0.20 httr_1.4.4 [127] AnnotationDbi_1.60.0 Rdpack_2.4 [129] colorspace_2.1-0 XML_3.99-0.13 [131] tensor_1.5 reticulate_1.28 [133] splines_4.2.2 statmod_1.5.0 [135] uwot_0.1.14 spatstat.utils_3.0-1 [137] scater_1.26.0 sp_1.6-0 [139] plotly_4.10.1 xtable_1.8-4 [141] jsonlite_1.8.4 nloptr_2.0.3 [143] R6_2.5.1 pillar_1.8.1 [145] htmltools_0.5.4 mime_0.12 [147] glue_1.6.2 fastmap_1.1.0 [149] minqa_1.2.5 BiocParallel_1.32.5 [151] BiocNeighbors_1.16.0 codetools_0.2-19 [153] utf8_1.2.3 lattice_0.20-45 [155] spatstat.sparse_3.0-0 numDeriv_2016.8-1.1 [157] pbkrtest_0.5.2 ggbeeswarm_0.7.1 [159] leiden_0.4.3 gtools_3.9.4 [161] survival_3.5-3 glmmTMB_1.1.5 [163] repr_1.1.6 munsell_0.5.0 [165] GetoptLong_1.0.5 GenomeInfoDbData_1.2.9 [167] iterators_1.0.14 variancePartition_1.28.0 [169] reshape2_1.4.4 gtable_0.3.1 [171] rbibutils_2.2.13
rc.integrated <- readRDS("../../CheWei/scRNA-seq/Integrated_Objects/rc.integrated_6S_GL2_lines_seu3_20240423.rds")
rc.integrated
An object of class Seurat 71741 features across 37615 samples within 3 assays Active assay: SCT (24847 features, 0 variable features) 2 other assays present: RNA, integrated 3 dimensional reductions calculated: pca, umap, umap_2D
table(rc.integrated$orig.ident)
sc_130 sc_131 sc_133 sc_134 sc_135 sc_137 6656 7658 8424 7837 4975 2065
rc.integrated$geno <- rc.integrated$orig.ident
rc.integrated$geno <- gsub("sc_130","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_134","WT",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_131","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_135","bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_133","pGL2_BRI1_GFP_bri1_T",rc.integrated$geno)
rc.integrated$geno <- gsub("sc_137","pGL2_BRI1_GFP_bri1_T",rc.integrated$geno)
rc.integrated$rep <- rc.integrated$orig.ident
rc.integrated$rep <- gsub("sc_130","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_134","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_131","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_135","2",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_133","1",rc.integrated$rep)
rc.integrated$rep <- gsub("sc_137","2",rc.integrated$rep)
table(rc.integrated$orig.ident, rc.integrated$geno)
bri1_T pGL2_BRI1_GFP_bri1_T WT
sc_130 0 0 6656
sc_131 7658 0 0
sc_133 0 8424 0
sc_134 0 0 7837
sc_135 4975 0 0
sc_137 0 2065 0
feature_names <- read_tsv("../data/features.tsv.gz", col_names = c("AGI", "Name", "Type")) %>%
select(-Type) %>%
distinct()
Rows: 32833 Columns: 3 ── Column specification ──────────────────────────────────────────────────────── Delimiter: "\t" chr (3): AGI, Name, Type ℹ Use `spec()` to retrieve the full column specification for this data. ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rc.integrated$geno <- factor(rc.integrated$geno, levels=c("WT", "bri1_T", "pGL2_BRI1_GFP_bri1_T"))
time_zonecell_typetime_zone_cell_typetime_zone_cell_subtypesrc.integrated$cell_type <- rc.integrated$celltype.anno.Li.crude
rc.integrated$time_zone <- rc.integrated$time.anno
rc.integrated$time_zone_cell_type <- rc.integrated$time.celltype.anno.Li.crude
table(rc.integrated$orig.ident, rc.integrated$cell_type)
Atrichoblast Columella Cortex Endodermis Lateral Root Cap Pericycle
sc_130 1468 893 982 520 1078 271
sc_131 1405 669 800 595 1787 749
sc_133 1421 605 1452 884 1684 332
sc_134 1195 588 924 639 2136 665
sc_135 1040 287 541 525 963 407
sc_137 351 282 293 236 421 134
Phloem Procambium Trichoblast Xylem
sc_130 57 190 1021 176
sc_131 119 331 944 259
sc_133 54 145 1669 178
sc_134 123 353 937 277
sc_135 106 211 707 188
sc_137 27 68 197 56
table(rc.integrated$time.anno)
Distal Columella Distal Lateral Root Cap Elongation
2943 2941 10498
Maturation Meristem Proximal Columella
7025 8801 359
Proximal Lateral Root Cap
5048
# Plot celltype annotation Li
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400D3", "#DCD0FF","#5AB953", "#BFEF45", "#008080", "#21B6A8", "#82B6FF", "#0000FF","#E6194B", "#DD77EC", "#9A6324", "#FFE119", "#FF9900", "#FFD4E3", "#9A6324", "#DDAA6F", "#EEEEEE")
rc.integrated$cell_type <- factor(rc.integrated$cell_type, levels = order[sort(match(unique(rc.integrated$cell_type),order))])
color <- palette[sort(match(unique(rc.integrated$cell_type),order))]
options(repr.plot.width=15, repr.plot.height=7)
(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'geno', pt.size = 0.75, ncol=3))
ggsave("../output/GL2/pGL2_BRI1-GFP_cell_types.pdf", width=13, height=6)
options(repr.plot.width=10, repr.plot.height=6)
FeaturePlot(rc.integrated, features="GL2-BRI1-GFP", split.by = "geno", order=T)
#ggsave("../output/GL2/pGL2_BRI1-GFP_expression.pdf", width=10, height=6)
Warning message in FeaturePlot(rc.integrated, features = "GL2-BRI1-GFP", split.by = "geno", : “All cells have the same value (0) of GL2-BRI1-GFP.”
options(repr.plot.width=10, repr.plot.height=6)
FeaturePlot(rc.integrated, features="AT1G79840", split.by = "geno", order=T)
#ggsave("../output/GL2/GL2-AT1G79840_expression.pdf", width=10, height=6)
options(repr.plot.width=30, repr.plot.height=7)
(BR_cell <- DimPlot(rc.integrated, reduction = "umap", group.by = "cell_type", cols = color, split.by = 'orig.ident', pt.size = 0.75, ncol=8))
DefaultAssay(rc.integrated) <- "SCT"
FeaturePlot(rc.integrated, features="GL2-BRI1-GFP", split.by = "orig.ident", order=T)
Warning message in FeaturePlot(rc.integrated, features = "GL2-BRI1-GFP", split.by = "orig.ident", : “All cells have the same value (0) of GL2-BRI1-GFP.” Warning message in FeaturePlot(rc.integrated, features = "GL2-BRI1-GFP", split.by = "orig.ident", : “All cells have the same value (0) of GL2-BRI1-GFP.”
saveRDS(rc.integrated, file = "/hpc/group/pbenfeylab/CheWei/scRNA-seq/Integrated_Objects/rc.integrated_6S_GL2_lines_seu4_annotated_20240429.rds")